##################################################
#
# PLOTS OF HARSH REPLICATION MODELS FOR
# OIL, ISLAM, and WOMEN PROJECT
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Created 27 February 2016
#
##################################################

##################################################
# R Housekeeping
##################################################

# Remove and empty workspace
rm(list=ls())

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Packages
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

library(foreign)
library(sandwich)
library(apsrtable)
library(tonymisc)  # includes summaryR()
library(MASS) # for polr() which runs our ologit

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Functions
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Code",sep=""))

source("statifyFN.R")
source("multipleOLSFN.R")
source("genModelsFN.R")
source("plotlineFN.R")
source("plotprepFN.R")
source("modelDefinitions.R")

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Data
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 


setwd(paste(mywd,"/Data",sep=""))

data <- read.dta("NewDVModelsXSforISQFINAL2015.dta")

data$oil_gas_pc_63 <- data$oil_gas_pc_63/1000

#################################################
#
# Lagged Replication: FLFP
# 
#################################################

out.harsh.rossalt <- multipleOLS(genModels(c("laborfemaleMOD_2000"), 
                                           p140alt),modelnames=names(p140base))
# Note change in modelnames parameter from p140alt to p140base
out.harsh.yupaulalt <- multipleOLS(genModels(c("laborfemaleMOD_2000"), p140earlyalt),
                                   modelnames=names(p140early))
# Note change in modelnames parameter from p140earlyalt to p140early


## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
# Plotting Models
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
vars.to.plot <- c("oil_gas_valuePOPred","logGDPcap_2000","age15_64",
                  "logGDPcap6372", "wdi_agedr_63","oil_gas_pc_63")

plot.harsh.rossalt <- do.call(rbind.data.frame,lapply(out.harsh.rossalt,plotprep,vars.to.plot))
plot.harsh.rossalt <- data.frame(plot.harsh.rossalt,c("Ross"))
colnames(plot.harsh.rossalt)[ncol(plot.harsh.rossalt)] <- c("Author")

plot.harsh.yupaulalt <- do.call(rbind.data.frame,lapply(out.harsh.yupaulalt,plotprep,vars.to.plot))
plot.harsh.yupaulalt <- data.frame(plot.harsh.yupaulalt,c("Ours"))
colnames(plot.harsh.yupaulalt)[ncol(plot.harsh.yupaulalt)] <- c("Author")

oilvars <- c("oil_gas_valuePOPred","oil_gas_pc_63")

plot.harsh.oipc <- rbind(plot.harsh.rossalt[plot.harsh.rossalt$Var%in%oilvars,],
                         plot.harsh.yupaulalt[plot.harsh.yupaulalt$Var%in%oilvars,])
plot.harsh.oipc$Offset <- .1
plot.harsh.oipc[plot.harsh.oipc$Author=="Ross",]$Offset <- -.1	
plot.harsh.oipc$Index <- 1:(nrow(plot.harsh.oipc)/2)
plot.harsh.oipc$X <- plot.harsh.oipc$Index + plot.harsh.oipc$Offset
plot.harsh.oipc$pch <- 21
plot.harsh.oipc[plot.harsh.oipc$Author=="Ross",]$pch <- 24
plot.harsh.oipc$col <- "black"
plot.harsh.oipc$bg <- "black"
plot.harsh.oipc$bg <- ifelse(abs(plot.harsh.oipc$Estimate/plot.harsh.oipc$SE)>=1.96,"black",
                              ifelse(abs(plot.harsh.oipc$Estimate/plot.harsh.oipc$SE)<=1.65,"white",
                                     "gray"))
incomevars <- c("logGDPcap_2000","logGDPcap6372")
plot.harsh.gdppc <-  rbind(plot.harsh.rossalt[plot.harsh.rossalt$Var%in%incomevars,],
                           plot.harsh.yupaulalt[plot.harsh.yupaulalt$Var%in%incomevars,])
plot.harsh.gdppc$Offset <- .1
plot.harsh.gdppc[plot.harsh.gdppc$Author=="Ross",]$Offset <- -.1	
plot.harsh.gdppc$Index <- 1:(nrow(plot.harsh.gdppc)/2)
plot.harsh.gdppc$X <- plot.harsh.gdppc$Index + plot.harsh.gdppc$Offset
plot.harsh.gdppc$pch <- 21
plot.harsh.gdppc[plot.harsh.gdppc$Author=="Ross",]$pch <- 24
plot.harsh.gdppc$col <- "black"
plot.harsh.gdppc$bg <- "black"
plot.harsh.gdppc$bg <- ifelse(abs(plot.harsh.gdppc$Estimate/plot.harsh.gdppc$SE)>=1.96,"black",
                            ifelse(abs(plot.harsh.gdppc$Estimate/plot.harsh.gdppc$SE)<=1.65,"white",
                                   "gray"))

agevars <- c("wdi_agedr_63","age15_64")
plot.harsh.age <-  rbind(plot.harsh.rossalt[plot.harsh.rossalt$Var%in%agevars,],
                         plot.harsh.yupaulalt[plot.harsh.yupaulalt$Var%in%agevars,])
plot.harsh.age$Offset <- .1
plot.harsh.age[plot.harsh.age$Author=="Ross",]$Offset <- -.1	
plot.harsh.age$Index <- 1:(nrow(plot.harsh.age)/2)
plot.harsh.age$X <- plot.harsh.age$Index + plot.harsh.age$Offset
plot.harsh.age$pch <- 21
plot.harsh.age[plot.harsh.age$Author=="Ross",]$pch <- 24
plot.harsh.age$col <- "black"
plot.harsh.age$bg <- ifelse(abs(plot.harsh.age$Estimate/plot.harsh.age$SE)>=1.96,"black",
                            ifelse(abs(plot.harsh.age$Estimate/plot.harsh.age$SE)<=1.65,"white",
                            "gray"))


#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing OIPC
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd("../Drafts/Charts")
#pdf(file="ISQFINALHarshFLFPSimple.pdf",height=5,width=7)
tiff(file="ISQFINALMainBodyFigure3A.tiff",
		height=5,
		width=7,
		units='in',
		res=300)
par(mfrow=c(3,1),oma=c(0,0,3,1),mar=c(3,3,1,1))
plot(1,
     pch="",
     ylim=c(1.15*min(plot.harsh.oipc$Y0),-1.15*max(plot.harsh.oipc$Y1)), 
     xlim=c(.5,max(plot.harsh.oipc$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of Oil Income Per Capita",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.harsh.oipc$Y0,
         y1=plot.harsh.oipc$Y1,
         x0=plot.harsh.oipc$X,
         x1=plot.harsh.oipc$X,
         col=plot.harsh.oipc$col,
         lwd=1,
         lty=1)
points(x=plot.harsh.oipc$X,
       y=plot.harsh.oipc$Estimate,
       col=plot.harsh.oipc$col,
       pch=plot.harsh.oipc$pch,
       bg=plot.harsh.oipc$bg,
       cex=1.5)
axis(side=1,at=plot.harsh.oipc[1:8,]$Index,
     labels=plot.harsh.oipc[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)


#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing Age
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

plot(2,
     pch="",
     ylim=c(1.15*min(plot.harsh.age$Y0),1.15*max(plot.harsh.age$Y1)), 
     xlim=c(.5,max(plot.harsh.age$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of Working Age Population",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.harsh.age$Y0,
         y1=plot.harsh.age$Y1,
         x0=plot.harsh.age$X,
         x1=plot.harsh.age$X,
         col=plot.harsh.age$col,
         lwd=1,
         lty=1)
points(x=plot.harsh.age$X,
       y=plot.harsh.age$Estimate,
       col=plot.harsh.age$col,
       pch=plot.harsh.age$pch,
       bg=plot.harsh.age$bg,
       cex=1.5)
axis(side=1,at=plot.harsh.age[1:8,]$Index,
     labels=plot.harsh.age[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing Age
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

plot(3,
     pch="",
     ylim=c(1.15*min(plot.harsh.gdppc$Y0),1.15*max(plot.harsh.gdppc$Y1)), 
     xlim=c(.5,max(plot.harsh.gdppc$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of GDP Per Capita",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.harsh.gdppc$Y0,
         y1=plot.harsh.gdppc$Y1,
         x0=plot.harsh.gdppc$X,
         x1=plot.harsh.gdppc$X,
         col=plot.harsh.gdppc$col,
         lwd=1,
         lty=1)
points(x=plot.harsh.gdppc$X,
       y=plot.harsh.gdppc$Estimate,
       col=plot.harsh.gdppc$col,
       pch=plot.harsh.gdppc$pch,
       bg=plot.harsh.gdppc$bg,
       cex=1.5)
axis(side=1,at=plot.harsh.gdppc[1:8,]$Index,
     labels=plot.harsh.gdppc[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)


# Title for whole thing
mtext("Ross's FLFP Models vs Models Employing Severe Lags", side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()


#################################################
#
# The Harsh Replication: Women in Parliament
# 
#################################################

# Here we use a 2002 DV and compare our models of
# women's representation in parliament
# (using 1963-1972 values of explanators) with
# Michael's (using 1993-2002 values of explanators).

out.harsh.ross.parliament <- multipleOLS(formulas=genModels(c("female_seats_2000"), 							p142base),
					modelnames=names(p142base))
out.harsh.yupaul.parliament <- multipleOLS(formulas=genModels(c("female_seats_2000"), 						p142early),
 					modelnames=names(p142early))

## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
# Plotting Models
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
vars.to.plot.parliament <- c("oil_gas_valuePOPred","logGDPcap_2000",
	"logGDPcap6372", "me_nafr_2000","oil_gas_pc_63")
				
plot.harsh.ross.parliament <- do.call(rbind.data.frame,
                                      lapply(out.harsh.ross.parliament,plotprep,vars.to.plot.parliament))
plot.harsh.ross.parliament <- data.frame(plot.harsh.ross.parliament,c("Ross"))
colnames(plot.harsh.ross.parliament)[ncol(plot.harsh.ross.parliament)] <- c("Author")

plot.harsh.yupaul.parliament <- do.call(rbind.data.frame,lapply(out.harsh.yupaul.parliament,plotprep,vars.to.plot.parliament))

plot.harsh.yupaul.parliament <- data.frame(plot.harsh.yupaul.parliament,c("Ours"))

colnames(plot.harsh.yupaul.parliament)[ncol(plot.harsh.yupaul.parliament)] <- c("Author")


plot.harsh.oipc.parliament <- rbind(plot.harsh.ross.parliament[c(3,6,9,12),],
		plot.harsh.yupaul.parliament[c(3,6,9,12),])
plot.harsh.oipc.parliament$Offset <- .1
plot.harsh.oipc.parliament[plot.harsh.oipc.parliament$Author=="Ross",]$Offset <- -.1	
plot.harsh.oipc.parliament$Index <- 1:4
plot.harsh.oipc.parliament$X <- plot.harsh.oipc.parliament$Index + plot.harsh.oipc.parliament$Offset
plot.harsh.oipc.parliament$pch <- 21
plot.harsh.oipc.parliament[plot.harsh.oipc.parliament$Author=="Ross",]$pch <- 24
plot.harsh.oipc.parliament $col <- "black"
plot.harsh.oipc.parliament $bg <- "black"
plot.harsh.oipc.parliament$bg <- ifelse(abs(plot.harsh.oipc.parliament$Estimate/plot.harsh.oipc.parliament$SE)>=1.96,"black",
                            ifelse(abs(plot.harsh.oipc.parliament$Estimate/plot.harsh.oipc.parliament$SE)<=1.65,"white",
                                   "gray"))		
						
plot.harsh.gdppc.parliament <- 	rbind(plot.harsh.ross.parliament[c(1,4,7,10),],
						plot.harsh.yupaul.parliament[c(1,4,7,10),])
plot.harsh.gdppc.parliament$Offset <- .1
plot.harsh.gdppc.parliament[plot.harsh.gdppc.parliament$Author=="Ross",]$Offset <- -.1	
plot.harsh.gdppc.parliament$Index <- 1:4
plot.harsh.gdppc.parliament$X <- plot.harsh.gdppc.parliament$Index + plot.harsh.gdppc.parliament$Offset
plot.harsh.gdppc.parliament$pch <- 21
plot.harsh.gdppc.parliament[plot.harsh.gdppc.parliament$Author=="Ross",]$pch <- 24
plot.harsh.gdppc.parliament$col <- "black"
plot.harsh.gdppc.parliament$bg <- "black"
plot.harsh.gdppc.parliament[plot.harsh.gdppc.parliament$Author=="Ross",]$bg <- "white"
plot.harsh.gdppc.parliament$bg <- ifelse(abs(plot.harsh.gdppc.parliament$Estimate/plot.harsh.gdppc.parliament$SE)>=1.96,"black",
                                        ifelse(abs(plot.harsh.gdppc.parliament$Estimate/plot.harsh.gdppc.parliament$SE)<=1.65,"white",
                                               "gray"))  	

						
plot.harsh.mena.parliament <- rbind(plot.harsh.ross.parliament[c(2,5,8,11),],
						plot.harsh.yupaul.parliament[c(2,5,8,11),])
plot.harsh.mena.parliament$Offset <- .1
plot.harsh.mena.parliament[plot.harsh.mena.parliament$Author=="Ross",]$Offset <- -.1	
plot.harsh.mena.parliament$Index <- 1:4
plot.harsh.mena.parliament$X <- plot.harsh.mena.parliament$Index + plot.harsh.mena.parliament$Offset
plot.harsh.mena.parliament$pch <- 21
plot.harsh.mena.parliament[plot.harsh.mena.parliament$Author=="Ross",]$pch <- 24
plot.harsh.mena.parliament$col <- "black"
plot.harsh.mena.parliament$bg <- "black"
plot.harsh.mena.parliament$bg <- ifelse(abs(plot.harsh.mena.parliament$Estimate/plot.harsh.mena.parliament$SE)>=1.96,"black",
                                         ifelse(abs(plot.harsh.mena.parliament$Estimate/plot.harsh.mena.parliament$SE)<=1.65,"white",
                                                "gray"))    

setwd("../Drafts/Charts")

# FOR PDF
# pdf(file="ISQFINALHarshWomenParliamentSimple.pdf",height=5,width=7)

# FOR TIFF (HI-RES)
tiff(file="ISQFINALMainBodyFigure3B.tiff",
		height=5,
		width=7,
		units='in',
		res=300)


par(mfrow=c(3,1),oma=c(0,0,3,1),mar=c(3,3,1,1))

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing OIPC
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

plot(1,
     pch="",
     ylim=c(1.15*min(plot.harsh.oipc.parliament$Y0),1.15*max(plot.harsh.oipc.parliament$Y1)), 
     xlim=c(.5,max(plot.harsh.oipc.parliament$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of Oil Income Per Capita",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.harsh.oipc.parliament$Y0,
		y1=plot.harsh.oipc.parliament$Y1,
		x0=plot.harsh.oipc.parliament$X,
		x1=plot.harsh.oipc.parliament$X,
		col=plot.harsh.oipc.parliament$col,
		lwd=1,
		lty=1)
points(x=plot.harsh.oipc.parliament$X,
		y=plot.harsh.oipc.parliament$Estimate,
		col=plot.harsh.oipc.parliament$col,
		pch=plot.harsh.oipc.parliament$pch,
		bg=plot.harsh.oipc.parliament$bg,
		cex=1.5)
axis(side=1,at=plot.harsh.oipc.parliament[1:8,]$Index,
		labels=plot.harsh.oipc.parliament[1:8,]$Model,
		cex.axis=1.25, tick=TRUE)
		

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing Age
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

plot(2,
     pch="",
     ylim=c(1.15*min(plot.harsh.mena.parliament$Y0),1.15*max(plot.harsh.mena.parliament $Y1)), 
     xlim=c(.5,max(plot.harsh.mena.parliament $Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of Middle East Regional Dummy",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0= plot.harsh.mena.parliament $Y0,
		y1= plot.harsh.mena.parliament $Y1,
		x0= plot.harsh.mena.parliament $X,
		x1= plot.harsh.mena.parliament $X,
		col= plot.harsh.mena.parliament $col,
		lwd=1,
		lty=1)
points(x= plot.harsh.mena.parliament $X,
		y= plot.harsh.mena.parliament $Estimate,
		col= plot.harsh.mena.parliament $col,
		pch= plot.harsh.mena.parliament $pch,
		bg= plot.harsh.mena.parliament $bg,
		cex=1.5)
axis(side=1,at= plot.harsh.mena.parliament[1:8,]$Index,
		labels= plot.harsh.mena.parliament[1:8,]$Model,
		cex.axis=1.25, tick=TRUE)
			
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing Age
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
				
plot(3,
     pch="",
     ylim=c(-1.15*min(plot.harsh.gdppc.parliament$Y0),1.15*max(plot.harsh.gdppc.parliament $Y1)), 
     xlim=c(.5,max(plot.harsh.gdppc.parliament $Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimates of GDP Per Capita",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0= plot.harsh.gdppc.parliament $Y0,
		y1= plot.harsh.gdppc.parliament $Y1,
		x0= plot.harsh.gdppc.parliament $X,
		x1= plot.harsh.gdppc.parliament $X,
		col= plot.harsh.gdppc.parliament $col,
		lwd=1,
		lty=1)
points(x= plot.harsh.gdppc.parliament $X,
		y= plot.harsh.gdppc.parliament $Estimate,
		col= plot.harsh.gdppc.parliament $col,
		pch= plot.harsh.gdppc.parliament $pch,
		bg= plot.harsh.gdppc.parliament $bg,
		cex=1.5)
axis(side=1,at= plot.harsh.gdppc.parliament[1:8,]$Index,
		labels= plot.harsh.gdppc.parliament[1:8,]$Model,
		cex.axis=1.25, tick=TRUE)
		

# Title for whole thing
mtext("Ross's Parliamentary Models vs Models Employing Severe Lags", side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()
